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Abstract 

Using Monte Carlo simulation, Van Duijneveldt and Lekkerkerker [Phys. Rev. Lett. 71, 4264 
(1993)] found gas-liquid-solid behaviour in a simple two-dimensional lattice model with two 
types of hard particles. The same model is studied here by means of numerical transfer matrix 
calculations, focusing on the finite size scaling of the gaps between the largest few eigenvalues. 
No evidence for a gas-liquid transition is found. We discuss the relation of the model with a 
solvable RSOS model of which the states obey the same exclusion rules. Finally, a detailed 
analysis of the relation with the dilute three-state Potts model strongly supports the tricritical 
point rather than a three-phase point. 



1 Introduction 

The phase behaviour of hard particles, in particular spheres, as a simple model of interacting 
particles, has received much attention. Computer simulations of monodisperse hard spheres show a 
first-order transition between a dilute disordered phase (fluid) and a dense ordered phase (solid) ||l|, 
^, ^. The continuous translational symmetry of the Hamiltonian remains intact in the fluid, but 
is broken to a discrete subgroup in the solid. Although a rigorous proof is lacking, this phase 
transition in the hard sphere model is now generally accepted. For bidisperse hard spheres the 
situation is more complicated. The existence of several solid phases has been established; see for 
example [Q and the references therein. The behaviour in the fluid phase, however, is not known. 
Using the Percus-Yevick closure of the Ornstein-Zernike equation, Lebowitz and Rowlinson |^ 
found miscibility in all proportions for all diameter ratios. More recently however, Biben and 
Hansen using the Rogers- Young closure, found a spinodal instability when the diameter ratio 
exceeds 5. Even so it might be that the fluid-fluid transition is pre-empted by the fluid-solid 
transition, so that the former does not actually occur. Thus it remains an open question whether 
bidisperse spheres can show a fluid-fluid phase separation. More generally one may ask if gas- 
liquid-solid behaviour can occur in binary mixtures with only hard-core repulsion. 

Motivated by this interest Van Duijneveldt and Lekkerkerker 0, |^ studied a two-dimensional 
binary hard-core lattice model. This model, introduced by Frenkel and Louis [|| , consists of large 
and small hard hexagons on a triangular lattice, see Figure]^. Every site can be empty or occupied 
by a large or small hexagon, and if it is occupied with a large hexagon all its direct neighbours 
must be empty. When the small particles are omitted, one regains the hard hexagon model [|o|, 
which has been solved exactly by Baxter |2j; it has a second-order ordering transition. Van 
Duijneveldt and Lekkerkerker studied the binary model by means of Monte Carlo simulation. 
They found three phases: dilute disordered (gas), dense disordered (liquid), and ordered (solid). 
Figure ^ shows this phase diagram, represented in terms of the fugacities zi and Z2 of the large 
and small hexagons, respectively. 

In this paper we study the same model by different methods. Our interest is in the quali- 
tative, rather than quantitative, aspects of the phase diagram. We do not address the general 
question whether gas-liquid-solid behaviour is possible in binary hard-core mixtures. The paper 



Figure 1: A typical configuration of large and small hexagons. 




Figure 2: Phase diagram in the 21-2:2 plane calculated by Van Duijneveldt and Lekkerkerker |^ 
from Monte Carlo simulations. The letters F, G, L and S indicate the fluid, gas, liquid and soUd 
phase, respectively. 



is organised as follows. First we briefly review the Monte Carlo approach of Van Duijneveldt and 
Lekkerkerker, and we give some exact results. Then we describe our numerical transfer matrix 
calculations. Next we discuss the relation of the model with an exactly solvable RSOS model and 
with the dilute three-state Potts model. Finally we propose an explanation for the discrepancy 
between our results and those of Van Duijneveldt and Lekkerkerker. 

2 Monte Carlo simulation and exact results 

Before we review the Monte Carlo method of Van Duijneveldt and Lekkerkerker ^ and discuss 
some exact results, we make the following notational conventions: the subscripts 1 and 2 refer 
to the large and small hexagons, respectively; the superscript refers to the pure hard hexagon 
model; the symbol N without subscript is the number of sites and is generally omitted as an 
argument of the thermodynamic quantities. 

We consider the semi-grand canonical partition function Z{Ni,Z2) of large hexagons, whose 
number TVi is fixed, and small hexagons, whose fugacity Z2 is fixed, on N lattice sites. We may view 
the small hexagons as causing an effective so-called depletion interaction between the large 
hexagons. The question is then if this attractive depletion interaction is strong enough to induce 
a fluid-fluid transition. The effective interaction can be expressed in the number of sites available 
for small hexagons, once the large hexagons have been placed on the lattice. Interestingly, the 
sites available for small hexagons are exactly the sites where an additional large hexagon could be 
inserted. Such sites are called free. It is easy to express the semi-grand canonical partition function 
Z{Ni, Z2) in terms of the canonical partition function Z^{Ni) of the hard hexagon model and the 
probability distribution p(A^f |iVi) for the number Nf of free lattice sites in the hard hexagon model: 

Z(^i,Z2) =Z°(iVi)Z'(A^i,Z2), 

where 

Z'(7Vi, Z2) = ^ p(iVf |iVi)(l + za)"^'-. 

After taking logarithms this gives the free energy: 

F{Ni,Z2) = F\Ni) + F'{Ni,Z2), (1) 

Van Duijneveldt and Lekkerkerker determine the probability distribution p from canonical Monte 
Carlo simulations of the hard hexagon model. To determine accurately the wings of the distribution 
an umbrella sampling technique is employed. They calculate F' from p, and for fixed Z2 fit a 
polynomial in pi :— Ni/N to this quantity. They obtain the free energy F from (^, using 
Baxter's exact result for and the fitted polynomial for F'. The fugacity zi of the large 

hexagons and the pressure P are calculated in the usual way from F. Finally phase equilibrium 
is determined by looking for phases with equal zi and P but different pi. As this calculation is 
carried out for fixed Z2, Z2 is also equal in the phases. The resulting phase diagram is shown in 
Figure ||. It has three branches: liquid-solid, gas-solid and gas-liquid. The branches meet at 
the three-phase point, at zi = 22.5 and Z2 — 1.89. (Van Duijneveldt and Lekkerkerker use the 
term "triple point", but as that suggests the coexistence of three phases where three first-order 
transitions meet we prefer to use the term "three-phase point".) The gas-liquid end-point is 
located at zi = 13.3 and Z2 = 1.36. 

Expanding Z' to first order in Z2 gives: 

Z'{Ni,Z2) = 1 + Z2{Ni)%^ + 0{Z2). (2) 

For a finite system we could have written 0(z|) instead of 0(22), but in the thermodynamic limit 
this is not valid at the phase transition of the hard hexagon model. Lekkerkerker (unpublished) 
found that the average pi : = {Ni/N)^_^ can be calculated exactly, as follows. Adding one hexagon 
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to a configuration of A'^i hexagons can be done in A^f ways. By doing tliis to all configurations of 
Ni hexagons each configuration of iVi + 1 hexagons is obtained exactly A^i + 1 times. Hence 



{Ni)%^Z°{m) = (TVi + 1)Z"{N, + 1), 
which in the thermodynamic limit yields 

- -. (3) 

This is an example of Widom's famous particle- insertion formula . In the Appendix we apply 
this exact result in the method of Van Duijneveldt and Lekkerkerker. In particular we show that 
the existence of a Van der Waals loop cannot be concluded from its presence in the first order 
approximant (|^). 

As the first derivatives of the thermodynamic functions with respect to Z2 are known in this way, 
we shall now attempt to calculate the locus of the phase transition in this order. The difference 
between the large and small hexagons is, that two small hexagons may occupy neighbouring sites, 
whereas two large ones may not. At small Z2 the density of small hexagons is low, so that they will 
generally occur isolated. Thus they cannot be distinguished from the large ones. For the grand 
canonical partition function this implies: 

Z{ZI,Z2) = Z"{zi+Z2)+0{Z2). (4) 

This suggests that the locus of the phase transition is given by 

Zl = - Z2 + 0(Z2), (5) 

where the superscript c refers to the critical point of the pure hard hexagon model. The particle 
densities follow also: 

Pl{zi, Z2) = ; p\{zi + Z2) + o{z2) 

Zl + Z2 

for the large hexagons, and similarly for the small ones. Combining these results yields the density 
of the large hexagons at the phase transition: 

Pi=(l-f)jp\ + o{z2). (6) 

Equations (||) and cannot be derived rigorously from (||) alone, but we conjecture that they 
are nevertheless valid. 



3 Transfer matrix approach 

Now we study the model through its row-to-row transfer matrix. For practical reasons, we work 
with sawtooth rows as shown in Figure ||. One advantage is that the high density ground state of 
the hexagons fits on the lattice (which has an even number of sites), whereas for straight rows it 
does so only when the system size is a multiple of three. Another advantage is that the transfer 
matrix can be built up by repeatedly adding one site, without increasing the total number of sites. 
Periodic boundary conditions are imposed on the rows. The number of "teeth" is denoted by W , 
so a row contains 2W sites and has length L = W^/2>. The largest few eigenvalues of the transfer 
matrix (in the zero-momentum sector) were calculated numerically for W = 2, . . . , b, using the 
power method. 

In the ordered regime there are in fact three coexisting ordered phases, corresponding to the 
three sub-lattices of the triangular lattice. They give rise to three eigenvectors of the transfer 
matrix, dominated by these ordered phases: one symmetric and two asymmetric for permutations 
among the ground states. The symmetric vector has the largest eigenvalue Aq. The asymmetric 
vectors have a complex conjugate pair of eigenvalues Am and A£^. In the region of the phase 
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Figure 3: The transfer matrix adds one row (shaded) to the system. 



diagram we are interested in, Aq, Am and AJ^, together with another real eigenvalue At, turn out 
to be the largest eigenvalues of the transfer matrix. The phase behaviour can be diagnosed from 
the behaviour of the gaps between the eigenvalues. Am : = log |Ao/Am| and At : = log |Ao/At|, 
as the system size L tends to infinity. 

The gap At is an inverse correlation length between density fluctuations. In the absence of 
a phase transition, the bulk (L = oo) value of this length is finite and the value for finite L 
approaches this bulk value when L tends to infinity. Hence At tends to a non-zero limit. At a 
critical point the bulk correlation length diverges and the value for finite L is proportional to L. 
As a consequence of scale invariance At decreases as At a first-order transition with a 

change in the density, however. At is not an inverse correlation length. The eigenvalues Ag and 
At are then asymptotically degenerate. Their gap At is related to the interfacial tension between 
the coexisting phases. More precisely. At decays as exp(— cri), where a is proportional to the 
interfacial tension [ p5[ . 

For the gap Am the situation is analogous. In the disordered regime, it is an inverse correlation 
length, here between fluctuations in the sub-lattice ordering. Thus the gap approaches a non-zero 
value as L grows. At a first-order transition between two disordered phases this correlation length 
is generally different in the two phases. Therefore the value of Am undergoes a sharp change 
through the transition, approaching a jump as the system size L increases. At a critical point the 
bulk correlation length diverges, so that Am decays as 1/L when L increases. In the ordered regime 
three phases coexist, and the eigenvalues Aq and Am (and A^) are asymptotically degenerate: Am 
decays exponentially with L. At a first-order transition between an ordered and a disordered phase 
by the same token Am vanishes exponentially with L. 

We shall now distinguish between two scenarios: (i) there are two phases (fluid and solid) as 
in Figure ^(a); (ii) there are three phases (gas, liquid and solid) as in Figure §(b). The gaps 
should behave as follows. At fixed 22, the gap Am decreases with increasing zi, whereas At has 
a minimum at the phase transition(s). For low Z2, see the lower dashed lines in Figures ^(a) and 
^(b), the scaled gaps I/Am and LAt will tend to a non-zero value when L — > 00 at the critical 
line. For high Z2, see the upper dashed lines, this is no longer the case: both scaled gaps tend to 
zero when L ^ 00 at the phase transition, which is now first-order. On the middle dashed line in 
Figure ^(b). Am changes rapidly at the gas-liquid transition. Furthermore At has two minima: 
at the gas-liquid transition and at the liquid-solid transition. When L — s- cx), the minimum of 
the scaled gap LAt tends to zero at the gas-liquid transition, but to a non-zero value at the 
liquid-solid transition. Thus the gas-liquid transition in Figure ^(b) can be recognised from the 
appearance of a sudden change in Am and a second minimum of At. 

For Z2 = 0.0, 0.1, . . . , 3.0 the scaled gaps LAm and LAt were plotted as function of zi 
for W = 2 . . , 5. Figures ||-^ show examples of this. We found no indication that At has 
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Figure 4: (a) Phase diagram with a fluid and sohd phase. The critical hne (fat) terminates at a 
tricritical point where the phase transition becomes first-order (double line), (b) Phase diagram 
with gas, liquid, and solid phases. The critical line (fat) meets the first-order transition (double 
line) at the three-phase point. 

two minima. One could argue that two minima might be fused to a single one for these relatively 
small systems; however, the sharpest and deepest minimum (at the gas-liquid transition) is clearly 
absent. This pleads against the three-phase scenario in favour of the two-phase scenario. We also 
saw no sudden change in Am. However, even if a gas-liquid transition were present, the signal in 
Am might be hard to detect. 

The three-phase scenario can be obtained by introducing an extra parameter into the model. 
Assign a weight k, to every lattice edge joining a small hexagon and an empty site. For k = 1 
one recovers the original model. For k = any contact between a small particle and an empty 
site is forbidden. In this limit the model either contains no small hexagons at all or is completely 
filled with them. The regime without small hexagons still exhibits the hard hexagon transition 
as long as 1 -|- Z2 is smaller than the partition sum per site of the hard hexagon model. Beyond 
this value the phase filled with small particles takes over. Thus the ordered and disordered hard 
hexagon phases meet with the pure small hexagon phase, where the phase transition between them 
terminates in a three-phase point. For k close to zero the model will still obey the three-phase 
scenario. Here At is indeed found to have two minima, see Figure ^. (The maxima in this figure 
at first sight seem to be crossings of eigenvalues, but a very close look reveals that they are in 
fact rounded.) This supports our interpretation of the absence of a second minimum in Ax as 
evidence against the three-phase scenario. 

The locus in the 21-2:2 plane of the phase transition can be estimated for example as the 
location of the minimum of Ax- For fixed 22 the value of 21 at which this gap takes its minimum 
was determined. The results iov W — h and W — 6 are plotted in Figure In order to obtain 
the locus in the pi~22 plane the density of large hexagons was computed using 

Pi = 21^ (-logAo) . 
021 

(It should be noted that for such small W this does not seem to be very accurate.) Figure |ll| 
shows the result. We observed that for fixed 22 the graphs of pi versus 21 for different system sizes 
pass approximately through one point. One could ask whether this is the critical point, as would 
be the case in a self-dual model. The locus of the intersection of the graphs for W = 5 and W = 6 
is shown in Figure Figures |l^ and ^ also show the phase diagrams given by Van Duijneveldt 
and Lekkerkerker |q|7 
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Figure 9: The scaled gaps LAt as a function of Zi on the hne Z2 = 1-3 in the model with extra 
parameter k — 0.6. The inset shows the deep minima in more detail. 




Figure 10: Locus in the ^1-2:2 plane of the minimum of the gap At for W = 5 (+) and W = 6 
(x) and phase diagram of Van Duijneveldt and Lekkerkerker (solid line). The asymptote (||) is 
also shown. 
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Figure 11: Locus in the pi-Z2 plane of phase transition calculated from W = 5 (+) and W = 6 
(x), locus of the intersection of the graphs ior W = 5 and W = 6 of pi versus zi (©), and phase 
diagram of Van Duijneveldt and Lekkerkerker (solid line). The asymptote (||) is also shown. 



First-order and second-order transitions are not easily distinguished from each other by the 
numerical data. In both cases At has a minimum, only the dependence on L of the depth of 
the minimum is different. For Z2 — 1.7, the graphs of the LAm pass approximately through one 
point, see Figure ^. The LAt have a minimum that increases slowly with L and may converge to 
a non-zero value, see Figure ^. This points to a second-order transition. For Z2 = 2.3, the graphs 
of LAm do not pass neatly through one point, see Figure |^. The minimum of LAt decreases 
with L and may vanish asymptotically, see Figure |8[ This points to a first-order transition. The 
behaviour of LAm and LAt changes gradually between Z2 = 1.7 and Z2 = 2.3. Thus the value of 
Z2 at the tricritical point is estimated roughly to lie between 1.7 and 2.3. 

By universality the limit values of LAm and LAt at the phase transition are 27rxM and 27ra;T 
respectively, with xm = 2/15 and xt — 4/5 on the hard hexagon critical line (c = 4/5), and 
xm = 2/21 and xt = 2/7 at the hard hexagon tricritical point (c — 6/7), see for instance [ p^ . 
On the critical line close to the critical point one expects to find the tricritical values for small 
system sizes, but the critical values for large sizes. The limits were also estimated from the graphs 
of LAm and LAt for Z2 = 0.0 (not shown) and Z2 — 1.7. For Z2 = 0.0 we found xm ~ 0.14 and 
XT ~ 0.80. This is in good agreement with the critical values xm = 2/15 and xt = 4/5. For 
Z2 = 1.7 we found xm ~ 0.13 and xt ~ 0.3. This agrees reasonably with the tricritical values 
Xm = 2/21 and xt = 2/7, which are expected for small system size near the tricritical point. 



4 Relation to an A2 RSOS model 

Some properties of the large-and-small hexagon model it has in common with an exactly solvable 
model. In order to make use of the exact solution we investigate if the two models are ever 
parametrically close. The sites of the large-and-small hexagon model can be in three states: 

(empty), 1 (large hexagon), or 2 (small hexagon). For neighbouring sites the combinations 1-1 

(2) 

and 1-2 are excluded. The same is true for the L = 7 case of the exactly solvable A2 restricted 
solid-on-solid model of Kuniba ||l^, |l^ . This is an interaction-round-a-face model on the square 
lattice. For a suitable choice of its spectral parameter, the condition on neighbouring sites extends 
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to one of the diagonals of the square face. The Boltzmann weight of the square face then factors 
into weights of the composing triangles: 



and these triangle weights are invariant under rotation: 

"'(: »)="f «)^"'(c «)^"'C °)^"'(^ I 

so that the model is isotropic on the triangular lattice. The model still has one parameter (the 
elliptic nome), but this solvable line stays away from our phase diagram. For example at the 
critical point the triangle weights are 

-(S o) ' ^ 

O' = 4-412, 
w{^^ Q I = 3.903, 
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2 2 



which is not of the form 
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Application of the numerical transfer matrix method from Section ^ to this critical model shows 
that it is in the tricritical three-state Potts universality class. 



5 Relation to the dilute three-state Potts model 

The large-and-small hexagon model is intimately related to the dilute three-state Potts model [|9| . 
Because this relation gives insight in the phase diagram we will consider it here in more detail. 
On every site j of a two-dimensional lattice with coordination number v lives a variable Sj that 
can take the values 0, 1, 2, 3. Of these the states Sj > take the role of local occupancy of one of 
the three sub-lattices of the hard hexagon model, and the state Sj = is neutral or vacant. The 
Hamiltonian of the dilute Potts model is 

E {Ss„s,+KSs^.oSs,,o)~lY,Ss„o, (7) 

<j,k> j 
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where the first sum is over nearest neighbour pairs of sites. In the parameter space {K,L,T) the 
model has a fine of tricritical points as well as a line of critical-end-points , see Figure As 
we will argue below, it is fairly clear where these come together, namely in the critical point of the 
four-state Potts model, K = 0, L = and T = T^, where all the four states are treated identically. 




Figure 12: A qualitative picture of the phase diagram of the dilute three-state Potts model. The 
dense coexistence region (back) and the dilute region (front) are separated by the three-state Potts 
critical surface (shaded) and the lower part of the first-order surface (not shaded). These surfaces 
meet at a line of three-state Potts tricritical points (left) and a line of three-state Potts critical- 
end-points (right). The upper part of the first-order surface (not shaded) separates a dilute and a 
dense disordered phase. It is bounded by a line of Ising critical points. The fat dot indicates the 
four-state Potts critical point. 

At T = there is a dilute phase with Sj — when vK + 2L > 0, while the three dense, 
or ordered phases associated with Sj — 1,2,3 coexist when vK + 2L < 0. These phases extend 
to non-zero temperatures so that a first-order surface separates the dilute region from the dense 
coexistence region. This first-order surface will not remain precisely at vK + 2L — for T > 0, but 
by symmetry it does include the T-axis, K = L ~ 0. At high temperature the coexistence region 



is bounded by a surface of three-state Potts critical points, shaded gray in Figure 12, where the 
line tension between the coexisting dense phases vanishes. This critical sheet must join with the 
first-order surface in a line of multicritical points, as they both form boundaries to the coexistence 
region. 

The nature of this multicritical line depends on the sign of K, as follows. Along the first-order 
sheet we can distinguish two line tensions, namely that between two different dense phases, and 
that between a dense and the dilute phase. When K <0 the interface between the dilute and the 
dense phases costs less energy than that between two of the dense phases. However on the critical 
surface the line tension between the dense phases vanishes. As a consequence all line tensions 
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vanish simultaneously where the critical and first-order sheets meet as K < 0. The separatrix 
between these two types of phase transition is thus a tricritical line. When K > the dense-dense 
interface costs less energy than the dense-dilute interface, so there remains a positive line tension 
between the dilute phase and the dense phases where the first-order sheet meets the critical surface, 
and the dense-dense interfacial tension vanishes. This results in a critical-end-point scenario: The 
three-state Potts critical sheet terminates where it hits the first-order sheet. The first-order sheet 
extends beyond this line, separating a disordered dense phase from the dilute phase. Obviously, 
at if = the two scenarios come together, and we conclude that the tricritical curve and the 
critical-end curve as well as the critical line terminating the dilute-disordered phase transition 
all meet in the four-state Potts critical point, marked as a dot in Figure |2|. This qualitative 
description of the phase diagram of (^, though not rigorous, is the simplest possible scenario, and 
has been corroborated by numerical studies . 

These considerations are of interest for the large-and-small hexagon model because that can be 
mapped onto a model sufficiently similar to the dilute Potts Hamiltonian (0) that the arguments 
can be carried over. We divide the triangular lattice into triangular blocks of three sites each, 
indicated in Figure ^(a). Each block then has three sites which we label 1, 2, and 3. We assign 
a spin variable Sj to each block, as follows. When the site a in block j is occupied by a large 
hexagon, the spin variable takes the value Sj — cr, while in all other cases Sj — 0. For convenience 
of notation we consider one block variable sq, in interaction with six neighbours Sj with 1 < |j | < 3, 
as shown in Figure ^(a). The blocks j with j > contain two sites neighbouring the site j of 
the central block, and the block — j sits in the opposite direction. To give an expression for the 
interaction we introduce the variables 



P^ = {Ss,,0 + Ss,.^) (1 - Ss_^,j) (l - (5,_,,fe) , (8) 

where i, j,k is a permutation of 1, 2, 3. Note that pi can only take the values and 1, and it signals 
if site i of the central block is free. The spin states 1, 2, and 3 have weight zi, but are excluded 
by some configurations of the neighbouring blocks by the factor 

[l-5s„A^-p,)]. (9) 

In other words the state sq = j is not allowed when pj = 0. The weight of the spin state sg = 
depends on the surrounding blocks and is given by the expression 

{1 + Z2)'"+P'+P' . (10) 

If this model would be precisely the dilute Potts model with Hamiltonian (|^) we could simply 
read off the value of K and its sign would conclusively decide between a tricritical point versus 
a three phase point. The interaction is of course much more complicated than that of the dilute 
Potts model, but the overall effect is that some combinations of unequal nearest neighbours are 
excluded or suppressed. As the state is treated altogether different from the states 1, 2, and 3, 
it is difficult to judge the sign of the effective couphng K in (^. 

However, this problem can be resolved because there is a model in the universality class and 
with the symmetry of the four-state Potts model which can be mapped to a very similar model. 
Consider a one-species lattice gas on the triangular lattice in which not only first neighbours but 
also second neighbours (at distance \/3) can not be occupied simultaneously. We will refer to this 
model as the big hexagon model. For large values of the fugacity z this model will be in an ordered 
phase in which one out of four sub-lattices is occupied preferentially. At low fugacity the symmetry 
between the sub-lattices is unbroken. The phase transition is known to be in the four-state Potts 
universality class from the symmetry of its Landau-Ginzburg- Wilson Hamiltonian ^ . We 
are not aware of studies giving the critical fugacity of this model, but we have seen numerically 
that it is about half the value of the hard hexagon model. 

The big hexagon model can be mapped exactly onto a Potts-like model very similar to the model 
above, as expressed in (^) and (|l^). Now we take blocks of four sites as shown in Figure p^(b), 
one in each sub-lattice. It is convenient to label the spins in each block by the numbers 0, 1, 2, 
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(a) (b) 

Figure 13: (a) The large-and-small hexagon model can be mapped onto a Potts- hke model by 
grouping the sites into blocks of three. The numbers indicate the labelling of blocks and of the 
sites within the blocks, (b) The big hexagon model can be mapped onto a Potts-like model by 
dividing the sites into blocks of four. The numbers indicate the labelling of the sites within the 
blocks. The blocks are numbered as in (a). 

3 as indicated. When the site j in a block is occupied, the block variable takes the value j. In 
addition, when none of the sites are occupied, the block variable is taken to be 0. Therefore the 
weight of the states j > is z and the weight of state will again depend on the states of the 
neighbouring blocks. We again consider a block variable sq interacting with its neighbours, which 
are labelled in the same way as in the previous case. We will use again variables pi defined by 
(^. The central site of the block is free if and only ii pi = p2 = P3 = 1- Some combinations 
of states of neighbouring blocks are excluded, described by precisely the same expression (||) as 
before. However, also some combinations of next-neighbouring blocks are excluded. For example 
site j of block — j and site k of block — fc in Figure |l^(b) are second neighbours, so the combination 
S-j = j and s_fe = fc is excluded. We introduce a variable 

Note that q can only take the values and 1; it signals if there are no pairs S-j = j and S-k = k. 
If So 7^ then s_j — j or s-k = fc is already excluded by the interaction between the neighbouring 
blocks and —j or — fc. Therefore the exclusion of the combination s-j — j and s_fc = fc can be 
taken into account by including a factor q in the weight of block in state 0. This weight is then 
given by 

q{l + z)P'P''P\ (11) 

In this way any exclusion between sites of next-neighbouring blocks is absorbed in the weight of 
state of the intervening block. 

This resulting model is strikingly similar to the Potts-like model above. The exclusion rules 
for pairs of neighbouring blocks are identical and when we choose Zi = z, the weight of the 
spin states 1, 2, and 3 is the same. In both models the weight of the state depends on the 
configuration of its six neighbours, via expression (^0|) and (pi]), respectively. When we further 
specify (1 + Z2)^ = (1 + z) the weights for sq — are equal in the case that pi — P2 — P3 
and q = I. In particular they are equal when the surrounding blocks are also in state 0, because 
then pi = P2 — P3 — I and q — 1- 
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It is the exclusion and suppression of configurations with unequal neighbours which determines 
an effective temperature T and coupling K in (^. The large-and-small hexagon model and the big 
hexagon model with the parameters as set above will have the same effective temperature T, as 
all configurations involving only spin states s > have the same weight between the two models. 
Only when a block has s = 0, while one or more of its neighbours have s > 0, the configurational 
weights between the two models can be different. In all such cases the weight in the big hexagon 
model is smaller than that in the large-and-small hexagon model, which is easy to see from direct 
comparison of the expressions (10) and (pi]). Therefore we can confidently claim that the effective 
coupling K is the greater in the big hexagon model, as configurations with unequal neighbours 
of which one s = are more strongly suppressed than in the large-and-small hexagon model. 
However, since the big hexagon model has the symmetry of the four-state Potts model, clearly 
its effective coupling K — Q. Therefore the effective K in the large-and-small hexagon model is 
necessarily negative, which, as argued above, results in a tricritical scenario. 



6 Discussion 

The results of our transfer matrix calculations provide evidence against the three-phase scenario 
of Figure ^(b) in favour of the two-phase scenario of Figure ^a). This contradicts the earlier 
findings of Van Duijneveldt and Lekkerkerker [Q, We propose the following explanation. Van 
Duijneveldt and Lekkerkerker effectively calculate the free energy difference between the binary 
mixture and the pure hard hexagons. They then look for phases of equal pressure and fugacities 
but different composition. They do not calculate the order parameter for the mixture. Their 
method has some drawbacks. Firstly, it cannot detect second-order transitions, because these 
do not involve a jump in the particle densities. Secondly, it uses a polynomial fit for the free 
energy difference, so that the total free energy still seems to possess the singularity of the pure 
hard hexagon model. Thirdly, whether P exhibits a Van der Waals loop or not may depend 
sensitively on p[Nf\Ni). Thus the locus of the liquid-solid branch in their phase diagram is a 
spurious consequence of the implicit assumption that the ordering transition remains at fixed pi 
for small values of Z2. Their qualitative conclusion that a gas-liquid transition is present relies on 
quantitative properties of the calculated phase diagram, viz. the locations of the various branches. 
Figure |l^ suggests that their gas-liquid and gas-solid branch together form the true fluid-solid 
line and that the critical point of their gas-liquid branch is in fact the tricritical point. This agrees 
well with the fact that Figures ^ and ^ show enhanced size dependence of the phase diagram near 
their gas-liquid critical point. However, this point is located at Z2 = 1.36 (and zi = 22.5), whereas 
we estimate roughly 1.7 < Z2 < 2.3 for the tricritical point. Being unable to present a satisfactory 
explanation for this discrepancy, we stress that our data do not signal a clearly determined locus 
of the tricritical point. It should also be noted that in our transfer matrix calculations only very 
small system sizes have been considered. Going to significantly larger systems might allow for 
more definitive quantitative statements, but this requires much greater computational resources. 

Other evidence comes from the relation with the dilute three-state Potts model. The large-and- 
small hexagon model can be mapped onto a Potts-like model. Another model, the big hexagon 
model, whose phase behaviour is known, can also be mapped onto a Potts-like model. A com- 
parison of the effective temperature and coupling constants between the large-and-small hexagon 
model on the one hand and the big hexagon model on the other hand indicates that the large- 
and-small hexagon model should follow the two-phase scenario. 
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Appendix 



It is instructive to follow the method of Van Duijneveldt and Lekkerkerker using (Q) and (||) instead 
of Monte Carlo results. Calculating the pressure from (p|) gives 



+ (ft-^)^2+o(z2). (12) 



Baxter p. 451] lists expansions around the critical point of several thermodynamic quantities 
of the pure hard hexagon model. Combining these expansions with (|) and (|l|) yields 

P = |p^ + ^^^||^sgn(pi-p5)|pi-p5r+0((pi-p5)^)| 

This suggests that for small non-zero values of Z2 the pressure P would exhibit a Van der Waals 
loop, so that the transition becomes first-order as soon as Z2 becomes non-zero. That this argument 
is not valid can be seen by considering for example 

f,{x) = {x-zf, 

which we view as a function of x, parametrically dependent on z. Expanding / to first order in z 
gives 

f,{x)^x'' -?,x^z + o{z) 

and for all non-zero values of z the function x^ — ix^z of x is decreasing between a; = and x — 2z. 
It is, however, a first order approximant of fz{x), which for all values of z is an increasing function 
of X. 
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